options("digits" =4)
coul <- c("#E0E0E0", "#FFFF99", "#CCFF99", "#00CC00", "#66FFFF","#0066CC", "#006633", "#FF99FF", "#994C00","#FF9999","#FF0000")
coul11 <- c(coul, "#FAE1C7", "#000000", "#FF8000", "#FAA4A4", "#d18975", "#2d543d", "#ff0055", "#5DA5DA", "#B276B2", "#F17CB0", "#009E73", "#F0E442", "#D55E00") #11 couleurs en ajouter en fct du nm de SM
my_palette = c("D0"= "#E0E0E0",
"D4"= "#ccFF99",
"D7"= "#0066CC",
"D10"= "#006633",
"D14"= "#FF99FF",
"D21"= "#994C00",
"D32"= "#FF9999",
"D90"= "#FF0000"
)
theme_CA <- function (base_size = 11, base_family = "", base_line_size = base_size/22,
base_rect_size = base_size/22)
{
theme_bw(base_size = base_size, base_family = base_family,
base_line_size = base_line_size, base_rect_size = base_rect_size) %+replace%
theme(panel.border = element_rect(fill = NA, colour = "black",
linewidth = rel(1)), panel.grid.major = element_blank(),
plot.title = element_text(size=10),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black", linewidth = rel(1)), legend.key = element_blank(),
strip.background = element_rect(fill = "white", colour = "black",
linewidth = rel(2)), complete = TRUE)
}
seu <- readRDS("/Users/victor/Documents/Master_2/Projet_long/Flux_metabo/results/seu__flux.RDS")
metabolite <- read.csv("/Users/victor/Documents/Master_2/Projet_long/Flux_metabo/data/Results_scFEA/mouse_balance_gut.csv", header = T, row.names = 1)
metabolite <- data.matrix(metabolite)
metabolite0 <- t(metabolite)
seu[['METABOLITES']] <- CreateAssayObject(counts = metabolite0)
DefaultAssay(seu) <- 'METABOLITES'
seu
## An object of class Seurat
## 29314 features across 30626 samples within 4 assays
## Active assay: METABOLITES (70 features, 0 variable features)
## 3 other assays present: RNA, SCT, FLUX
## 5 dimensional reductions calculated: pca, umap, pca.flux, umap.metabo, umap.flux
Profils de concentrations
day <- seu@meta.data$orig.ident
metabolite <- as_tibble(metabolite)
metabolite$day <- day
metabolite <- relocate(metabolite, day, .before = AMP)
dt_c0 <- as.matrix(seu@assays$METABOLITES@counts)
yyy <- seu@meta.data$orig.ident
levels(as.factor(seu@meta.data$orig.ident))
## [1] "D0" "D4" "D7" "D10" "D14" "D21" "D32"
coul_rev <- rev(coul)
for (jj in 1:nrow(dt_c0)) {
xxx <- dt_c0[jj, ]
final_df <- data.frame(var = paste('X', 1:length(xxx), sep=''), # le même que plus haut
meta = xxx,
day = yyy)
order <- rev(c("D0", "D3", "D4", "D5", "D6", "D7", "D10", "D14", "D21", "D32", "D90"))
final_df <- transform(final_df, day = factor(day, levels = order))
title <- rownames(dt_c0)[jj] # récupère le noms du module
aa <- ggplot(final_df, aes(x = meta, y = day, fill = day)) +
geom_density_ridges() +
theme_ridges() +
theme(legend.position = 'none') +
ggtitle(title) +
theme(plot.title = element_text(hjust = 0.5)) +
scale_fill_manual(values = coul_rev) +
xlab("Concentrations") +
geom_vline(xintercept = median(final_df$meta))
plot(aa)
}
## Picking joint bandwidth of 1.21e-05

## Picking joint bandwidth of 0.00141

## Picking joint bandwidth of 0.00019

## Picking joint bandwidth of 0.000509

## Picking joint bandwidth of 6.22e-05

## Picking joint bandwidth of 1.46e-07

## Picking joint bandwidth of 5.92e-07

## Picking joint bandwidth of 5.74e-06

## Picking joint bandwidth of 0.00111

## Picking joint bandwidth of 1.47e-06

## Picking joint bandwidth of 2.82e-07

## Picking joint bandwidth of 0.000479

## Picking joint bandwidth of 9.4e-08

## Picking joint bandwidth of 3.7e-06

## Picking joint bandwidth of 0.000644

## Picking joint bandwidth of 0.000598

## Picking joint bandwidth of 3.11e-08

## Picking joint bandwidth of 9.42e-06

## Picking joint bandwidth of 3.55e-05

## Picking joint bandwidth of 0.000146

## Picking joint bandwidth of 0.000848

## Picking joint bandwidth of 5.91e-06

## Picking joint bandwidth of 1.07e-05

## Picking joint bandwidth of 1.53e-05

## Picking joint bandwidth of 0.00156

## Picking joint bandwidth of 3.28e-06

## Picking joint bandwidth of 0.0066

## Picking joint bandwidth of 5.83e-07

## Picking joint bandwidth of 0.00427

## Picking joint bandwidth of 4.94e-06

## Picking joint bandwidth of 4.86e-05

## Picking joint bandwidth of 2.18e-09

## Picking joint bandwidth of 0.000131

## Picking joint bandwidth of 9.73e-06

## Picking joint bandwidth of 5.42e-06

## Picking joint bandwidth of 0.000631

## Picking joint bandwidth of 1.88e-06

## Picking joint bandwidth of 3.74e-05

## Picking joint bandwidth of 4.53e-06

## Picking joint bandwidth of 1.59e-06

## Picking joint bandwidth of 2.16e-09

## Picking joint bandwidth of 4.83e-06

## Picking joint bandwidth of 1.02e-08

## Picking joint bandwidth of 0.00121

## Picking joint bandwidth of 3.62e-05

## Picking joint bandwidth of 8.67e-06

## Picking joint bandwidth of 0.000136

## Picking joint bandwidth of 0.000112

## Picking joint bandwidth of 6.58e-08

## Picking joint bandwidth of 9.84e-06

## Picking joint bandwidth of 2.55e-06

## Picking joint bandwidth of 0.00013

## Picking joint bandwidth of 1.56e-06

## Picking joint bandwidth of 6.86e-07

## Picking joint bandwidth of 0.000152

## Picking joint bandwidth of 0.00165

## Picking joint bandwidth of 7.17e-05

## Picking joint bandwidth of 1.88e-06

## Picking joint bandwidth of 0.00506

## Picking joint bandwidth of 0.000129

## Picking joint bandwidth of 0.000598

## Picking joint bandwidth of 1.56e-08

## Picking joint bandwidth of 2.89e-06

## Picking joint bandwidth of 3.23e-09

## Picking joint bandwidth of 6.1e-07

## Picking joint bandwidth of 2.11e-07

## Picking joint bandwidth of 9e-05

## Picking joint bandwidth of 9e-05

## Picking joint bandwidth of 1.06e-06

## Picking joint bandwidth of 0.000259

Cinétique
means <- metabolite %>%
group_by(day) %>%
summarise(across(everything(), mean, na.rm= TRUE)) %>%
t()
## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `across(everything(), mean, na.rm = TRUE)`.
## ℹ In group 1: `day = D0`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))
means <- cbind(rownames(means), means)
colnames(means) <- as.character(means[1,])
colnames(means)[1] <- "metabolites"
means <- means[-1,]
means <- as_tibble(means)
means
## # A tibble: 70 × 8
## metabolites D0 D4 D7 D10 D14 D21 D32
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 AMP -0.001827 -0.001888 "-0.00… "-0.… -0.0… -0.0… -0.0…
## 2 Pyruvate 0.02118 0.10085 "0.054… "0.0… 0.02… 0.03… 0.02…
## 3 Acetyl.CoA -0.01212 -0.02067 "-0.01… "-0.… -0.0… -0.0… -0.0…
## 4 Glutamate -0.002500 -0.003629 " 0.00… "-0.… -0.0… -0.0… -0.0…
## 5 X2OG 0.007885 0.010517 "0.008… "0.0… 0.00… 0.00… 0.00…
## 6 Oxaloacetate -0.0005159 -0.0005223 "-0.00… "-0.… -0.0… -0.0… -0.0…
## 7 Glycine 0.005076 0.005048 "0.005… "0.0… 0.00… 0.00… 0.00…
## 8 Succinate 0.002201 0.002276 "0.002… "0.0… 0.00… 0.00… 0.00…
## 9 UDP.N.acetylglucosamine -0.005280 -0.001113 " 0.02… " 0.… -0.0… -0.0… -0.0…
## 10 lysine 0.01197 0.01194 "0.011… "0.0… 0.01… 0.01… 0.01…
## # ℹ 60 more rows
maxi <- c()
for (i in 1:70) {
maxi <- c(maxi, max(as.numeric(means[i, 2:8])))
}
mini <- c()
for (i in 1:70) {
mini <- c(mini, min(as.numeric(means[i, 2:8])))
}
delta <- maxi - mini
means$conc_min <- mini
means$conc_max <- maxi
means$conc_delta <- delta
means
## # A tibble: 70 × 11
## metabolites D0 D4 D7 D10 D14 D21 D32 conc_min conc_max
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 AMP -0.0… -0.0… "-0.… "-0.… -0.0… -0.0… -0.0… -1.93e-3 -1.82e-3
## 2 Pyruvate 0.02… 0.10… "0.0… "0.0… 0.02… 0.03… 0.02… 2.12e-2 1.01e-1
## 3 Acetyl.CoA -0.0… -0.0… "-0.… "-0.… -0.0… -0.0… -0.0… -2.07e-2 -1.21e-2
## 4 Glutamate -0.0… -0.0… " 0.… "-0.… -0.0… -0.0… -0.0… -4.55e-3 2.28e-3
## 5 X2OG 0.00… 0.01… "0.0… "0.0… 0.00… 0.00… 0.00… 7.89e-3 1.05e-2
## 6 Oxaloacetate -0.0… -0.0… "-0.… "-0.… -0.0… -0.0… -0.0… -5.22e-4 -5.16e-4
## 7 Glycine 0.00… 0.00… "0.0… "0.0… 0.00… 0.00… 0.00… 5.05e-3 5.08e-3
## 8 Succinate 0.00… 0.00… "0.0… "0.0… 0.00… 0.00… 0.00… 2.20e-3 2.28e-3
## 9 UDP.N.acetylgluc… -0.0… -0.0… " 0.… " 0.… -0.0… -0.0… -0.0… -5.28e-3 2.74e-2
## 10 lysine 0.01… 0.01… "0.0… "0.0… 0.01… 0.01… 0.01… 1.19e-2 1.20e-2
## # ℹ 60 more rows
## # ℹ 1 more variable: conc_delta <dbl>
means <- arrange(means, desc(conc_delta))
ggplot(means, aes(conc_delta)) + geom_dotplot() + scale_x_log10()
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

means_2 <- filter(means, conc_delta > 1e-3) ## Reste 32 métabolites
Profils de concentrations finaux
T_F <- rownames(dt_c0) %in% means_2$metabolites
dt <- dt_c0[T_F == T,]
coul_rev <- rev(coul)
for (jj in 1:nrow(dt)) {
xxx <- dt[jj, ]
final_df <- data.frame(var = paste('X', 1:length(xxx), sep=''), # le même que plus haut
meta = xxx,
day = yyy)
order <- rev(c("D0", "D3", "D4", "D5", "D6", "D7", "D10", "D14", "D21", "D32"))
final_df <- transform(final_df, day = factor(day, levels = order))
title <- rownames(dt)[jj] # récupère le noms du module
aa <- ggplot(final_df, aes(x = meta, y = day, fill = day)) +
geom_density_ridges() +
theme_ridges() +
theme(legend.position = 'none') +
ggtitle(title) +
theme(plot.title = element_text(hjust = 0.5)) +
scale_fill_manual(values = coul_rev) +
xlab("Concentrations") +
geom_vline(xintercept = median(final_df$meta))
plot(aa)
}
## Picking joint bandwidth of 0.00141

## Picking joint bandwidth of 0.00019

## Picking joint bandwidth of 0.000509

## Picking joint bandwidth of 6.22e-05

## Picking joint bandwidth of 0.00111

## Picking joint bandwidth of 0.000479

## Picking joint bandwidth of 0.000644

## Picking joint bandwidth of 0.000598

## Picking joint bandwidth of 3.55e-05

## Picking joint bandwidth of 0.000146

## Picking joint bandwidth of 0.000848

## Picking joint bandwidth of 0.00156

## Picking joint bandwidth of 0.0066

## Picking joint bandwidth of 0.00427

## Picking joint bandwidth of 4.86e-05

## Picking joint bandwidth of 0.000131

## Picking joint bandwidth of 0.000631

## Picking joint bandwidth of 3.74e-05

## Picking joint bandwidth of 0.00121

## Picking joint bandwidth of 3.62e-05

## Picking joint bandwidth of 0.000136

## Picking joint bandwidth of 0.000112

## Picking joint bandwidth of 0.00013

## Picking joint bandwidth of 0.000152

## Picking joint bandwidth of 0.00165

## Picking joint bandwidth of 7.17e-05

## Picking joint bandwidth of 0.00506

## Picking joint bandwidth of 0.000129

## Picking joint bandwidth of 0.000598

## Picking joint bandwidth of 9e-05

## Picking joint bandwidth of 9e-05

## Picking joint bandwidth of 0.000259

Clustering
seu <- FindVariableFeatures(seu, selection.method = "vst", verbose = F)
## Warning in eval(predvars, data, env): NaNs produced
## Warning in hvf.info$variance.expected[not.const] <- 10^fit$fitted: number of
## items to replace is not a multiple of replacement length
seu
## An object of class Seurat
## 29314 features across 30626 samples within 4 assays
## Active assay: METABOLITES (70 features, 70 variable features)
## 3 other assays present: RNA, SCT, FLUX
## 5 dimensional reductions calculated: pca, umap, pca.flux, umap.metabo, umap.flux
seu <- ScaleData(seu, features = rownames(seu), assay = 'METABOLITES', verbose = F)
seu <- RunPCA(seu, features = VariableFeatures(object = seu), reduction.name = 'pca.metabo', verbose = F)
## Warning in irlba(A = t(x = object), nv = npcs, ...): You're computing too large
## a percentage of total singular values, use a standard svd instead.
## Warning: Cannot add objects with duplicate keys (offending key: PC_), setting
## key to 'pca.metabo_'
DimPlot(seu, dims = 1:2,reduction = "pca.metabo", group.by = "orig.ident", label = T)

ElbowPlot(seu, ndims = 30, reduction = 'pca.metabo') # 9 dim

Clustree and clusterization

seu <- RunUMAP(seu, dims = 1:9, assay = 'METABOLITES', reduction = 'pca.metabo', reduction.name = "umap.metabo", verbose = F)
## Warning: Cannot add objects with duplicate keys (offending key: UMAP_) setting
## key to original value 'umap.metabo_'
seu <- FindNeighbors(seu, dims = 1:9, reduction = 'pca.metabo', verbose = F, assay = 'METABOLITES')
seu <- FindClusters(seu , graph.name = 'METABOLITES_snn', verbose = F)
DimPlot(seu, reduction = "umap.metabo", label = T, label.box = T, group.by = "METABOLITES_snn_res.0.3")

DimPlot(seu, reduction = "umap.metabo", label = T, label.box = T, group.by = "orig.ident")

DimPlot(seu, reduction = "umap", label = T, label.box = T, group.by = "METABOLITES_snn_res.0.3")

Barplot
cell <- cbind(Cells(seu), seu@meta.data$orig.ident)
colnames(cell) <- c("barcode", "dataset")
day <- seu@meta.data$orig.ident
cell <- as_tibble(cell)
cell$day <- day
cell$rna_cluster <- seu@meta.data$SCT_snn_res.0.8
cell$flux_cluster <- seu@meta.data$FLUX_snn_res.0.5
cell$metabo_cluster <- seu@meta.data$METABOLITES_snn_res.0.3
unique(cell$flux_cluster)
## [1] 12 16 4 11 5 14 10 17 6 3 8 13 7 2 1 15 9 0
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
unique(cell$rna_cluster)
## [1] 6 0 14 11 4 9 15 1 8 19 17 18 7 16 3 12 2 5 10 13
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
unique(cell$metabo_cluster)
## [1] 3 9 6 12 4 7 2 11 0 10 1 8 5
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12
ggplot(cell, aes(x = day, fill = rna_cluster)) +
geom_bar() +
labs(title = "Repartition des cluster RNA en fonction des jours",
y = "RNA clusters") +
theme_CA() +
scale_fill_manual(values = coul11) +
theme(axis.title.x = element_blank())

ggplot(cell, aes(x = rna_cluster, fill = day)) +
geom_bar() +
labs(title = "Repartition des jours dans les clusters RNA",
y = "Jour") +
theme_CA() +
scale_fill_manual(values = my_palette) +
theme(axis.title.x = element_blank())

ggplot(cell, aes(x = day, fill = flux_cluster)) +
geom_bar() +
labs(title = "Repartition des cluster de flux en fonction des jours",
y = "Flux clusters") +
theme_CA() +
scale_fill_manual(values = coul11) +
theme(axis.title.x = element_blank())

ggplot(cell, aes(x = flux_cluster, fill = day)) +
geom_bar() +
labs(title = "Repartition des jours dans les clusters de flux",
y = "Jour") +
theme_CA() +
scale_fill_manual(values = my_palette) +
theme(axis.title.x = element_blank())

ggplot(cell, aes(x = day, fill = metabo_cluster)) +
geom_bar() +
labs(title = "Repartition des cluster de flux en fonction des jours",
y = "Flux clusters") +
theme_CA() +
# scale_fill_manual(values = coul11) +
theme(axis.title.x = element_blank())

ggplot(cell, aes(x = metabo_cluster, fill = day)) +
geom_bar() +
labs(title = "Repartition des jours dans les clusters de métabolite",
y = "Jour") +
theme_CA() +
scale_fill_manual(values = my_palette) +
theme(axis.title.x = element_blank())
